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ABSTRACT 

Whether volatiles can be entrapped in a background matrix composing planetary en- 
velopes and be dragged via convection to the surface is a key question in understanding 
atmospheric fluxes, cycles and composition. In this paper we consider super-Earths with an 
extensive water mantle (i.e. water planets), and the possibility of entrapment of methane 



in their extensive water ice envelopes. We adopt the theory developed by van der Waals fc 



Platteeuw (1959) for modelling solid solutions, often used for modelling clathrate hydrates. 



and modify it in order to estimate the thermodynamic stability field of a new phase, called 
methane filled ice Ih. We find that in comparison to water ice VII the filled ice Ih structure 
may be stable not only at the high pressures but also at the high temperatures expected at 
the core-water mantle transition boundary of water planets. 



1. INTRODUCTION 

The discoveries and characterization of planetary systems orbiting other stars has en- 
tered an exciting period when we are starting to gain access to observing the atmospheres of 
planets that are essentially solid in nature - high-density rocky or icy planets of 1 to 10 Earth 
masses. These planets, called collectively super-Earths, have been discovered in relatively 



large quantities, though only a handful have measured radii and masses so far ( Carter et al. 



2012, and references therein). The mean densities derived for these exoplanets reveal a range 
of possible bulk compositions, ranging from rocky with high iron content (e.g., Kepler-lOb, 



Batalha et al. 2011) to mini-Neptunes with high H/He fraction (e.g., Kepler-lld,e, Lissauer 



et al. 2011). One of these planets orbiting a nearby M-dwarf star, GJ1214b (Charbonneau 



et al. 2009), has been accessible to spectroscopic studies of its atmosphere with inferences 



to its composition (e.g.. Bean et al. 2011 Berta et al. 2012, and references therein). 
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At intermediate densities many of the super-Earths should represent a new type of 
planet, unknown in our Solar System, composed of a rocky core and a water envelope that 
exceeds 10% by mass. We will call them water planets; this paper studies the transport of 
volatiles inside them. Possible examples of water planets could be Kepler-llb, 18b, or 20b. 



Massive water planets were introduced conceptually by Kuchner (2003) and Leger et al. 



(2004). Fu et al. (2010) developed models of the interior dynamics of the water layers of 



water planets and concluded that materials released from the silicate-iron core should reach 
the surface, despite the high pressures at the core-mantle boundary and phase transitions. 
Materials would traverse the water mantle, composed of high pressure ice phases, and reach 
the surface on timescales of 0.1 to 100 Myr. These are convective overturn timescales: what 
are the actual materials (e.g., gases) that could be transported to the surface and affect 
atmospheric observables needs to be determined by a detailed study. Here we attempt a 
first step in this direction by considering entrapment of volatiles under the extremely high 
pressures inside water super-Earths. 

The entrapment in water ice of volatile gases, which are mostly hydrophobic species, 
occurs in structures known as clathrate hydrates: the guest molecules occupy cages of wa- 
ter molecules stabilized predominantly by repulsive interactions between them. Typical 
clathrate hydrates include gases such as methane, carbon dioxide, oxygen, and nitrogen. 
Volatile transport by clathrates has been studied in the context of Titan and the icy satel- 



lites in our solar system ( [Tobie et al. 2006 Halevy & Stewart 2008 Schubert et al. 2010 



Lunine et al. 2010 Sohl et al. 2010). Because clathrates are of practical terrestrial, as well 



as of astrophysical importance, much effort has been devoted to measuring their properties 



in the laboratory (see, e.g. Fortes & Choukroun 2010, and the references therein). Most 
experiments are carried out at pressures up to an order of 100 MPa, which is high enough to 
provide useful information for modeling Titan-sized bodies. 

Experiments at higher pressure show that methane clathrates undergo a transition to 



a more compact form called filled ice at pressures of around 2GPa (Loveday et al. 2001b) 



that can survive at even higher pressures above 86 GPa (Hirai et al. 2006). For comparison 



in a typical water super-Earth Fu et al. (2010) estimate that the pressures at the bound- 



ary between the ice mantle and the silicate core will be of the order of 100 GPa while the 
temperature will be of the order of 10^ K. 



The basic theory of clathrates that was developed by van der Waals & Platteeuw ( 1959 ) 



was later applied by Lunine & Stevenson ( 1985 ) to situations of astrophysical interest. Below 



we suggest how this theory may be extended to higher pressures and temperatures. We use 
this theory to estimate the stability regime of filled ice, and discuss the implications of this 
for volatile transport in super-Earths. 



Clathrates are crystals, whose lattice structure forms cells that act as cages for foreign 
molecules. The empty cage-like structure is usually thermodynamically unstable, but the 
captured foreign molecules (i.e. guest molecules) help stabilize the clathrate crystal (van 



der Waals & Platteeuw 1959). Clathrate hydrate is essentially water ice. In this case a 



framework of groups of four coordinated water molecules creates a cage-like lattice. The 
hydrogen-bonded water molecules are slightly distorted, however, from the tetrahedral angle 
of ordinary ice. There are two basic low pressure (< 1 GPa) forms of clathrate hydrates, 
referred to as structure I (SI) and structure II (SII). Each structure is composed of two types 
of cages, one small and one large. There is an additional hexagonal structure (SH) whose 
thermodynamic stability regime is roughly in the 1 — 2 GPa range (Hirai et al. 2001). 



Filled water-ice structures also have the capability of trapping volatiles within their hy- 
drogen bonded network, except that instead of the clathrate cages the volatiles are entrapped 
in channels within the water ice ( Loveday et al.|2001b ). All cage clathrates have cages whose 
diameters are much larger than the channels which connect them, while the channels occu- 
pied by the guest molecules in filled ice are of constant and smaller diameter along their 
length (see Figure [I]). Their densities also differ significantly: for methane clathrate SI the 
water-methane ratio is 5.75:1, while for methane filled ice it is 2:1. The increased guest-host 
and guest-guest intermolecular interaction may help explain the increased stability of the 



filled ice against high pressure compression and decomposition (Machida et al. 2007). Dia- 



mond anvil experiments, using methane as the guest molecule, suggest such structures may 



survive pressures beyond 86 GPa (Hirai et al. 2006). 



The ability of clathrate structures to incorporate relatively large amounts of volatile 
species under conditions where the pure volatiles are highly unstable in the condensed phase 
is a primary characteristic that makes them so important to solar system and exoplanet 
studies. Enclathrating volatiles in icy planetesimals results in their inclusion in planetary 
icy mantles. As the icy mantle grows its inner layers become pressurized, and at pressures 
higher than about 2 GPa even SH clathrate hydrates are destabilized. Experiments show 
that the rate of pressure increase controls the resulting water phase: if the pressure is raised 
in steps of 0.2 — 0.5 GPa every 1 — 2hr, the methane clathrate will decompose to ice VII 
and pure methane ice. However, if the sample of methane clathrate hydrate is kept at more 
than 3 GPa for 12 hr the molecules have enough time to rearrange themselves into a filled ice 



structure that can sustain very high pressures without decomposing (Loveday et al. 2001a). 



Thus it appears that the filled ice phase will be an important component of a planetary ice 
mantle, especially inside water super-Earths. 

In this paper we begin by revisiting the thermodynamic stability theory of clathrates 
and extend it to higher pressures and temperatures. We consider only methane as the guest 



molecule here. In section 3 we consider the case of methane filled ice. We conclude with 
numerical results and brief application to super-Earths. 





Fig. 1. — Comparison between methane cage clathrate and methane fiUed ice Ih. The structure of fihed 
ice Ih is viewed down the c-axis, perpendicular to the widened channels formed in this phase. In water 
ice Ih every hexagon is hydrogen bonded to its neighbouring hexagons, along the c-axis, in an alternate 
manner. In methane filled ice Ih three adjacent oxygens of a particular hexagon (plus signs) will bond to 
one neighbouring hexagon and the other three (minus signs) will bond to the other neighbouring hexagon. 
This difference results in a widening of the channels perpendicular to the c-axis in filled ice-Ih, therefore 
allowing for the accomodation of methane. The black balls in the filled ice Ih structure (right panel) denote 



the entrapped methane molecules. (Filled ice after Loveday et al. (2001b) with permission) 



THERMODYNAMIC STABILITY FIELD 



The basic theory of clathrate hydrates views their stability in terms of two phases: The 
ordinary water phase (e.g. ice Ih, liquid water, etc. - the a phase) and the empty clathrate 
hydrate structure (the /3 phase) -|- a gas of guest molecules. The chemical potentials of 



the two phases are equal on the thermodynamic equihbrium boundary between them. For 



clathrates this may be written as (van der Waals & Platteeuw 1959): 



f^Q = ^^Q + kT^Uilnh 




VKi 



IkC 



Ki 



1 + Ej fJCJ^ ^^^ 

Here the subindex Q refers to the molecule that makes up the host molecular network, which 
in our case is H2O. The second term on the RHS of Eq.Q represents the contribution of 
the guest molecules to the clathrate hydrate chemical potential. T is the temperature, k is 
Boltzmann's constant, z/j is the ratio between the number of i type cages to water molecules 
per cubic unit crystal, and finally, yxi is the probability that a guest molecule of type K 
occupies a clathrate cage of type i. This last function is given in Eq.dll) in terms of the 
volatile fugacity {fx) and its Langmuir constant {Cxi)- We will give further information on 



the Langmuir constant below and refer the reader to van der Waals & Platteeuw ( 1959 ) for 
a more detailed derivation. 

Since most high pressure experimental information is for methane-filled water-ices we 
restrict ourselves to the case of a single species of guest molecule model, and omit the 
summation over K. The equilibrium equations may therefore be written as: 
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kT ~ kT 
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+ ^z/iln(l -ycH4,i) 
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(4) 



Recalling that both the fugacity and the Langmuir constant are functions of temperature 
and pressure, we differentiate and combine the last set of equations to get: 
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(5) 



We introduce the thermodynamic relations: 

9 ( l^H20\ _ Hh^o 
kT^ 
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Id 



kT 
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(7) 



Where HH20 and V^jO represent the enthalpy and volume per water molecule, respectively. 
Inserting these thermodynamic relations into Eq.dsl) and solving for dP/dT gives: 
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Here we have omitted the index CH4 for convenience. For a homogeneous substance 
(i.e. composed only of water molecules) the chemical potential equals the Gibbs free energy 
per particle, G, which, combined with the assumption of constant temperature gives: 



dn = dG = VdP = kTd\nP 



(9) 



Here all extensive parameters are per particle. For the final equality on the RHS we have 
used the equation of state for an ideal gas, so that the last relation between pressure and 
Gibbs free energy is only valid for this case. When the gas is not ideal, we can retain the 
functional form if we replace the pressure with the fugacity, /, which acts as an effective 
pressure function to correct for the effect of intermolecular interactions and which obeys the 



following relation (Smith & Van Ness 1975): 



dfx = dG = VdP = kTdlnf 



From the last relation we have: 
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Where V is the volume per methane molecule. For high pressures we can solve Eqs.(Q and 
d2| for the case of an SI methane clathrate hydrate numerically. We find that fCi ~ 10^ ^ 1. 



A similar result was found by Lunine & Stevenson (1985). Inserting this numerical result 



together with Eq.(ll) into Eq.rt8| yields: 






(12) 



is: 



The definition of the Langmuir constant according to van der Waals & Platteeuw ( 1959 ) 

h 

kTCK[T) 
where hxi is the single cell canonical partition function for a K type guest molecule in an 
i type cage. C,k is the quantum number density function for a K type molecule in an ideal 



gas and is independent of pressure. The cell partition function depends on pressure both via 
the cell dimension and through the form of the guest-host potential, so: 
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Here again we have omitted the index K. We are left with estimating the derivative of the 
single cell partition function, for which we give the following quasi-classical form: 
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Here we have divided the Hamiltonian of the guest molecule into its separate kinetic and 
potential contributions. As a first order approximation we assume that the kinetic degrees of 
freedom of the entrapped molecule are unaffected by its inclusion in the water network. This 
is a common approximation in clathrate modeling, and experiment shows that in solidified 



form methane molecules rotate freely as in an ideal gas (Hazen et al. 1980). In this ap- 



proximation the kinetic degrees of freedom will contribute a function that is only a function 
of temperature, and the logarithm will cancel upon differentiation with respect to pressure, 
thus yielding a simplified form: 
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In clathrate hydrates the cages entrapping the volatiles are assumed spherical (e.g. Sloan 



1998 McKoy & Sinanoglu 1963) so that the integration is over a spherical cage. In filled 



ices the cages are actually cylindrical channels within the water ice structure (Loveday et al. 



2001b) so that we may write: 
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Where we assume the cylindrical water-made channel has radius a and since intermolecular 
potentials fall rapidly with increasing distance we limit the integration along the z coordinate 
to the finite values zi and Z2. 

Both the limits of integration and the potential interaction of the methane molecule with 
its surroundings, Wi, may depend on the pressure. Raman spectra of methane filled water 
ice shows increases in the attractive potential between methane and its water network host 
with pressure increases ( Machida et al.||2007 ). This was also suggested by Hirai et al. (2006). 
The rearrangement of the water network, from cage clathrate to the filled ice structure, 
reduces molecular distances by ~ 0.5 x 10^^ cm. It is also suggested, from intermolecular 



distance considerations, that the guest-host Lennard- Jones potential interaction estimated 
for clathrates is enhanced by weak hydrogen bonds between guest and host in the filled ice 
structure ( Loveday et al.||2001b ). Since the guest-host potential energy changes considerably 
with pressure, a good first approximation will be to consider only the change of Wi with 
pressure, therefore allowing us to insert the derivative with respect to pressure into the 
integrand, giving: 
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Finally, taking a spatial average for the partial derivative appearing in the numerator, 
we have: 
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Inserting this last relation into Eq.(12) yields: 






(20) 



Since the pressure varies by many orders of magnitude and the temperature does not, and 

since pressure has a dominant effect in determining the correction for the system non-ideality, 

we assume that: 

' df 



— - > 
dT dT 



We may further approximate the following 

d 



kTj: 



fi t:= In C. 
dT 



kTj: 



Vi 



kTj: 

i 



^ 1 
dT 



1 
T 



1 
kf 



^3^g-W^./fcT 



T 



{W,) - k 



(21) 



Where in the last step we have again averaged the intermolecular interaction of the methane 
molecule spatially over its surroundings. Inserting the last formula into the relation for 
dP/dT gives a Clausius-Clapeyron type equation of the form: 



dP_^ (<o - H^H.o) + E. ^. [t (W^) - k] 
'dT 



(22) 



<0 - VS^o + E. ^^{m)T) - E. ^^V 

This modified Clausius-Clapeyron type equation is a generalization of the equation given in 



Lunine & Stevenson ( 1985 ) in order to explain the behaviour of clathrates. We will use this 



formalism to obtain the thermodynamic field of stability for the filled water ice as well. For 
this purpose we need to evaluate the different terms appearing in this equation, which we 
do in the next section. 



APPLICATION TO HIGH PRESSURE 

3.1. Clathrate Hydrate 



Starting with clathrate hydrates, let us examine the numerator of Eq.(22). At low 



pressures and temperatures the a phase is water-ice Ih whose enthalpy is taken to be equal 



to the enthalpy of the empty hydrate (/3) phase (Lunine & Stevenson 1985). The potential 
of interaction Wi is attractive (negative) and of a Lennard- Jones type, so generally the 
numerator is negative. In the denominator, at low pressures, the gaseous volatile volume 
{V) is the dominant factor appearing with a minus sign. Therefore the derivative dP/dT 
is positive. The potential of interaction hardly changes with increasing pressure in this low 
pressure regime. 

As we increase the temperature, the pressure increases till we reach the melting point 
for water-ice Ih and the a phase now represents liquid water. Due to the enthalpy of fusion, 
the enthalpy difference appearing in the numerator is no longer negligible and causes a sharp 
increase in the absolute value of the numerator. This is manifested as a sharp increase 
in the derivative dP/dT. Every increase in temperature is now accompanied by a steeper 
increase in pressure. At room temperature at 5 MPa, methane gas already deviates from 
ideality enough so that we need to consider a second virial correction; at 10 MPa a third 
virial correction is required, and so on ( Hirschf elder et al.|1966). This means that the volume 



a methane molecule occupies in the gas decreases with pressure. At high enough pressure 
the volatile gas contracts enough so that the empty clathrate volume equals the hquid water 
volume plus the compressed methane volume and the derivative dP/dT diverges. 

Any further increase in pressure brings about a situation where the volume occupied by a 
clathrate water molecule is larger than the sum of the volumes occupied by a water molecule 
in the liquid phase and a methane molecule in the gas phase weighted by the hydration 
number (z/). The result is that the derivative dP/dT becomes negative and the high pressure 
stability limit of the clathrate hydrate is attained and the clathrate dissociates. This general 
type of clathrate behavior is shown in fig. |2| which we have derived by numerically solving 
the set of Eqs.g and g. 

Now let us consider what happens for the case of filled water-ice. This ice is formed in 
the laboratory at room temperature and at a pressure of ~ 2 Gpa. The structure is an ice Ih 
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water network which is distorted in such a way that the interconnecting channels widen to 
accommodate the methane molecules within the hydrogen bonded network. This distortion 
from the usual structure of water ice Ih is manifested by a change in the tetrahedral angles 



(Loveday et al. 2001b). The /3 phase will now describe the filled ice structure while the a 



phase represent either water ice VII or fluid water, appropriately. 

To obtain the thermodynamic stability field for filled water ice we shall require the 
temperature and pressure dependencies for both the different constituents' volumes and for 
the attractive potential between a methane molecule and its surroundings. In addition, the 
enthalpy difference between filled ice and fluid water ought be estimated for the case of 
stability with respect to fluid water. 



3.2. High Pressure Equations of State 



We start with the general relation: 

dV 



XdT — ndP 



(23) 



Where x ^"^^ '^ ^ire the volumetric thermal expansivity and isothermal compressibility re- 
spectively. We assume the bulk modulus, B, is linearly dependent on pressure and may be 
written as: 

B = - = Bo + BoP (24) 

K 



Combining Eqs.(24) and (23) yields after integration: 

-I /Bo 



V{T,P) = V{To,Po 



Bq + BqP 
Bo + BoPo 



exp\ I xiT,P)dT\ 



(25) 



If we keep the temperature constant at Tq we end up with the Birch-Murnaghan equation of 
state. By setting the reference temperature (Tq) to room temperature, we can assign to each 
solid constituent a proper value for its bulk modulus (i?o) and its pressure gradient {Bq) by 
using published room temperature hydrostatic compression experiments. We use the data 



published in Hemley et al. (1987) for water ice up to 128 GPa, for assigning a bulk modulus 



for ice VII. For assigning a bulk modulus for solid methane we use the data published in 



Hazen et al. (1980) for phase I and in Umemoto et al. (2002) for phases A and B up to 



37 GPa. Solid methane transforms at high pressures to a hexagonal close packed structure 
( Bini fc Pratesi|[l997 ), this phase was difficult to account for as we found no volumetric data 



for it, rather we extrapolated from the phase B data in Umemoto et al. (2002) to higher 
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pressures. For the filled ice structure we used room temperature, experimentally deduced, 



unit cell volumes up to 42 GPa by Hirai et al. ( 2003 ) 



For the thermal expansivity we have adopted the approach of Fei et al. (1993), who 



determined the thermal expansivity of water ice VII by fitting their volumetric experimental 



data to an equation of state of the form given above (Eq,25). These authors assumed the 
following dependence for x oii the pressure: 



x(r,p) = x(T,Po 



Bo + BoPo 



-V 



(26) 



Where xiT,Po) is taken to be a linear function of the temperature and rj is well fitted 
with a numerical value of 0.9. For the thermal expansivity of solid methane we adopt the 



formalism of Eq.(26) for its dependency on pressure and set the expansivity value at Pq to be 
^ K~^, according to experimental data given by Heberlein & Adams (1970). The thermal 
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expansivity of filled ice is not known. As filled ice is a combination of a hydrogen bonded- 
network and methane molecules between which there is van der Waals type attraction, we 
suggest its thermal expansivity to be intermediate between that for ice VII and for solid 
methane. 

For the equation of state of fluid water at high pressure (up to 100 GPa) we adopt the 



formalism derived using molecular dynamic simulations by Belonoshko & Saxena (1991). 



Although this formalism is inherently dependent on the type of model used for the water 
molecules' intermolecular potential it does show good agreement with recent experimental 



data for high pressure water fluid density (Goncharov et al. 2009). 



We shall now turn to evaluate the energy of interaction of a methane molecule with its 
surroundings in the water filled ice structure. 



3.3. The Interaction Potential W^ 



In order to build the thermodynamic stability field for methane filled water ice, we 
need to approximate how the potential of interaction, of the methane molecule with its sur- 
roundings, depends on the temperature and the pressure. As was shown by [Raghavendra 



& Arunan (2008) a hydrogen bond may form between the water electron poor hydrogen 



and the center of the methane tetrahedral face which is electron rich. It was further shown 
by these authors that the bond energy (Efyond) ^ot such an interaction is —6.3 kJ mol^^. 
For comparison a methane-methane van der Waals potential well is about —1.2 kJ mol^^ 



(Hirschfelder et al. 1966). It was also suggested by Loveday et al. (2001b) that weak hydro 
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gen bonds are formed in the filled water-ice structure, between the water network and the 
dissolved methane molecules. 

Let us consider a simple approximate model where with increasing pressure more of 
the tetrahedral faces, per methane molecule, create such hydrogen bonds with the water 
network. An increase in temperature will have the opposite effect. That the number of 
hydrogen bonds per molecule increases with pressure and decreases with temperature is 



known for water structures (e.g. Kalinichev & Bass 1994 Pattanayak et al. 2011). We may 



therefore write for the spatially averaged potential of interaction of a methane molecule with 
its surroundings: 

{W,) = n{T,P)E,,nci (27) 

Where n{T, P) is the number of hydrogen bonds between a methane molecule and the water 
network, at a given pressure and temperature, and is bounded between and 4, where in 
the upper limit all of the methane four tetrahedral faces are hydrogen bonded to the water 
network. 

We normalize n to give it a probability interpretation (i.e. What is the probability a 
bond will form at a given pressure-temperature point). We further assume a division into a 
temperature dependent and pressure dependent probabilities, thus: 

n{T,P) = n,{T)n2{P) (28) 

For the temperature dependent probability we assume a Boltzmann type form, of: 

To obtain the form for the pressure dependent probability (^2) we use the fact that a for- 
mation of a hydrogen bond is accompanied by a substantial penetration of the hydrogen 
bonding molecules into each other. For the case of the hydrogen bonding between water and 
methane the combined van der Waals radius of 2.90 x 10~^ cm, before the bonding, reduces 
to 2.47 X 10~* cm, after bonding (Raghavendra fc Arunan||2008). This interpenetration may 



therefore account for a 40% reduction in the crystal volume. This is considerable enough so 
that we may relate this interpenetration (to first approximation) to the solid compressibility, 
K. At low pressures only a few hydrogen bonds are formed and the solid is easily compressed, 
as many bonds are still ready to be formed. As the pressure increases more bonds form, per 
molecule, and it becomes more difficult to further compress the solid. Hypothetically, at a 
high enough pressure, all bonds per molecule are already formed and it is no longer possible 
to further compress the solid via the route of hydrogen bonding molecular interpenetration. 
Adopting this model, to first approximation, we can write for the molecular volume, V, the 
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following: 

V = UhbVhB + rinHBVnHB (30) 

Where uhb and n-nHB are the probabilities a methane molecule is fully hydrogen bonded 
to the water network or not hydrogen bonded at all respectively, and Vhb and VnUB are 
the molecular volumes associated with these two molecular situations respectively. If the 
number of hydrogen bonds indeed determines, to first approximation, the crystal volume, 
we may say: 

''- V\dP)^ UHBVHB+n^nHBVnHBK dP ^ ^"^^ dP ) ^^'^ 



Since uhb + n^uB = 1, an integration of Eq.(31) yields: 



1 — exp ( — /„ KdP ) 
n2{P) = UHBiP) = -^^^^ ^ (32) 



Vr,. 



Inserting Eqs.(32) and (29) into Eq.(27) gives for the spatially averaged potential energy, of 



a methane molecule with its surroundings, the following form: 



(^^) - 1 ^ .-4L.I/.T (1 - ^-^'n E.on. (33) 



1 + Q-i\Ei,o^d\/kT 

Its derivative with respect to pressure will therefore be: 

dWi\ \ _ iHpEbond _KoP 



dP J^/ 1 + e-4|-E.o„dl/fcT 



e-^°^ (34) 



Now that we have the approximate temperature and pressure dependencies for the terms 



appearing in Eq.(22) we may integrate it numerically to obtain the stability regime for the 
filled ice structure. 



3.4. Numerical Results and Discussion 



The numerical integration of Eq. ( 22 ) is shown in fig. pi As we have mentioned above, we 
consider it reasonable to assume the thermal expansivity of the filled ice to be intermediate, 
between that for solid methane (a van der Waals solid) and water ice VII. From the laboratory 
data (see subsection 3.2) we know the thermal expansivity for solid methane is an order of 



magnitude larger than that for ice VII. We integrate Eq.(22) assuming thermal expansivity 
for the filled ice two times, three times and five times larger than the experimental value for 
ice VII. As a limiting case we also solve assuming filled ice has a thermal expansivity equal 
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to that of water ice VII. The four curves are given in the figure. It is appropriate to note 



here that Sloan ( 1998 ) gives for a SI clathrate hydrate (also a combination of methane and 



water) a thermal expansivity some five times larger than the thermal expansivity for ice Ih. 



More recent experiments confine the latter ratio to be between two to four (Hester et al. 



2007). 



In addition, the parameter u, the hydration number, is important to Eq.(22), as defined 



earlier. Generally speaking, this parameter represents the ability to include volatiles in the 



water network. As the ratio between water and methane in the filled ice is 2 : 1 (Loveday 



et al. 2001b) then filled ice modeling requires a value of 1/2 for u. Given that, we find the 



packing efficiency for filled ice to be greater than for the case of separation to pure solid 
constituents (i.e. water ice VII and pure solid methane). The contribution of the potential 
of interaction (in comparison to packing efficiency considerations) is found to be small and 
restricted to the lower pressure regime of the stability field. 



Integrating Eq.(22) and following the dissociation curve for filled-ice, with respect to 



ice VII, a point of intersection with the melting curve for water ice VII is reached. Such an 
intersection is commonly referred to as a quadruple point. Up to the quadruple point the 
enthalpy difference between the a and /3 phases is relatively small, as both phases are solid. 
Continuing the integration beyond the quadruple point the a phase will now represent fiuid 
water. The enthalpy difference between the a and /3 phases therefore increases, and becomes 
dominated by the enthalpy of fusion of filled ice. Unfortunately, the enthalpy of fusion of 
methane filled-ice Ih is experimentally undetermined and we estimate its value using the 
enthalpy of fusion for pure solid water at high pressure. 



Goncharov et al. (2009) deduced experimentally the melting curves for water ice VII and 



for super-ionic water. By further determining experimentally the volume difference between 
solid and melt, along the melting curve, they were able to derive the enthalpy of fusion for 



pure water. Their reported error is approximately 50%. As explained in Goncharov et al. 



(2009) the enthalpy of fusion increases with pressure since the increase in pressure means 



the melting transition is between a molecular solid and a fluid whose molecules become 
ever more dissociated and ionized. At a pressure of about 47 GPa a branching occurs in 
the melting curve for pure water, due to the introduction of the super-ionic phase of water 
( Goncharov et al^|2005 ) . The result is a reduction in the melting curve gradient and a sharp 



decrease of the enthalpy change upon melting (Goncharov et al. 2009). This behaviour is 



clearly seen in flgj4l where the dashed curve (red) reproduces the results from Goncharov 



et al. (2009) and the dashed-dotted curve (blue) is a hypothetical enthalpy of fusion, for 



which we have arbitrarily reduced the point of branching in the melting curve to 27 GPa. 
The reason for the introduction of this hypothetical behaviour is that the enthalpy of fusion 
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we adopt is from experiments on a homogeneous water system and therefore it is only an 
approximation for our methane filled-ice system. Molecular dynamic simulations byllitaka &; 



Ebisuzaki (2003) demonstrate that filled-ice exhibits behaviours similar to those appearing 



in homogenous water systems, but at somewhat lower pressures. A branching of the melting 
curve is therefore probable in our system of methane mixed with water as well, though, 
the branching may happen at a lower pressure then the 47 GPa of the pure water system. 
Wishing to test to what extent do our results depend on the point of branching we solve for 
the hypothetical enthalpy of fusion as well. 



In figJS^we solve Eq.(22) where we assume for filled-ice a thermal expansivity twice that 



for water ice VII and solve for each of the enthalpy of fusion scenarios depicted in figjlj In 



the same figure we also show the solution of Eq. ( 22 ) with the experimental enthalpy of fusion 



from Goncharov et al. (2009) but vary it globally by 50%, which is its experimental error. 



From the figure we see that most of the variation in the filled-ice Ih dissociation curve due 
to the changes examined in the enthalpy of fusion occur below 10 GPa. 

Combined with thermodynamic profiles for the interior of water exo-Earths, the ther- 
modynamic stability field can help us decide whether filled ice structures can form at the 
core-mantle boundary and help convect methane and other hydrocarbons towards the sur- 
face. 



APPLICATION TO SUPER-EARTHS 



For the interior structure of water planets with masses up to 10 M^, Fu et al. (2010) 



consider a silicate-metal core, surrounded by a water-ice mantle. They find sohd-state con- 
vection to prevail between the core-ice boundary and the conductive crust. Examining 
different surface temperatures and heat fluxes the authors find it possible for a liquid water 
ocean to exist beneath an ice-Ih solid surface layer. Different compounds expelled from the 
silicate-metal core may be incorporated in the ice matrix and convect outward. For the 
2M® planet Fu et al. (2010) find the sihcate core-ice boundary pressure to range from ap- 



proximately 50 to 90 GPa, for water mass fractions between 25% and 50% respectively. The 
expected temperatures are between 700 and 1000 K. 



Our extension of the theory of van der Waals & Platteeuw ( 1959 ) suggests that the filled 



ice structure may be stable at such pressures and temperatures. In such a case any CH4 
released from the core could be trapped in the lower part of the ice mantle. Experimentally 



Loveday et al. (2001b) inferred a 2 : 1 water-methane ratio for the filled ice phase. As the 



mantle convects, the CH4 would be carried to lower pressures where the filled ice will undergo 
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a transition to a clathrate hydrate. Using the statistical model given by [van der Waals fc 



Platteeuw (1959) we can compute the probability, y^i, that a structure I clathrate hydrate 
cage will be occupied. This is shown in fig.|6| It is clear that the occupancy is very close to 
full occupancy (probability of unity) for which there are 5.75 water molecules per methane 
molecule. For a structure H methane hydrate the water-methane ratio ranges from 4.25:1 to 
3.40:1, depending on the number of methane molecules occupying its large cage ( Koh||200"2 ). 



Thus as the pressure decreases in an upwelling, excess CH4 will be forced out of the water 
network. This may lead to the formation of a local methane reservoir. Such reservoirs may 
naturally occur on the transition from filled ice to structure H clathrate hydrate at about 
2 GPa, and at the transition from structure H clathrate hydrate to structure I clathrate 
hydrate at about 1 GPa. 

In A. Levi et al.(2013, in preparation) we show that the introduction of filled ice as a 
major constituent in an icy mantle has a large effect on the mantle thermodynamic profile. 
The probable higher thermal expansivity of filled ice compared with that for water ice VII 
results in a more moderate adiabatic thermal profile. The higher temperatures in the icy 
mantle, compared with those for a pure water mantle, creates a physical route through which 
super-Earths (objects less massive than Uranus and Neptune) may develop lower mantles 
in the super-ionic and reticulating phases. Phases so far related with the interior of bodies 
whose mass is equivalent to that of the icy giants of our solar system. 

The unique characteristics of methane clathrate hydrate, namely its low thermal con- 
ductivity and the topology of its melting curve, yield water planets with thin crusts (<1 km) 
and a tendency to form a layer closely confined to the local melting condition beneath it. In 
that respect the geology and surface-atmosphere coupling in water super-Earths are quite 
different then simply assuming a water planet is a scaled up version of an icy satellite, such 
as Titan. Icy satellites tend to form thick crusts, in the order of 100 km, and stabilize in a 
stagnant lid regime. This also means that chemical and geocycles in super-Earths are not 
simply resolved using analogies on icy moons. 

The formation of a layer, in which conditions are close to melting, underneath the 
crust acts as a low viscosity layer that may decouple the convection cell and the crust, 
and represent a channel for relatively fast horizontal flow from up-welling to down-welling 



(Crowley & O'Connell 2012). Such a fast horizontal flow may preferentially drag methane 
reservoirs to areas of down-welling where it may again become incorporated into the water 
matrix, preventing it from reaching the surface and the atmosphere. The time scale for 
transport of methane to the atmosphere is therefore broken down to several stages: the filled 
ice mantle overturn, the rate of transport in the low viscosity layer underlying the crust and 
the rate of collapse of the thin crust followed by exposure of fresh material. 
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Other substances trapped in filled ice structures may be dragged preferentially along 
with the convecting ice. If they reach the top of the convection cell and meet with the 
bottom of a subsurface ocean, some compounds may favour clathrate-hydrate formation. 
This will deplete those substances relative to others that are more easily transported in the 
liquid water and released to the atmosphere via cryovolcanism and surface ice tectonics. The 
volatile composition and abundance may also effect the ocean freezing temperature. This 
will change the thickness and temperature of the solid crust. These issues are developed in 
more detail in A. Levi et al.(2013, in preparation). 

In this work we have estimated the thermodynamic stability field for methane filled ice. 
From this field of stability it is clear that the blanketing effect of a thick H/He atmosphere, 
resulting in high internal temperatures, will most likely destabilize any filled ice (see figjs] 
for the P-T conditions at the base of the H/He atmosphere surrounding Uranus and Nep- 
tune). Water super-Earths lacking a substantial H/He atmosphere are more favourable for 
the formation of filled ice, which in turn will have a large influence on them. Additional 
experimental and theoretical work is needed to further constrain the properties of methane 
filled ice in order to better understand this influence. 

We wish to thank the anonymous referee for his constructive comments. 
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Fig. 2. — This figure describes the thermodynamic stabihty field for CH4 clathrate hydrate structure I. 
The dotted fine is water ice Ih melting curve, the thick curve profile (coloured blue in the online version) is 
the dissociation pressure curve for CH4 clathrate hydrate structure I, which clearly follows the available data 
points (circles). The other melting curves are for water ice III, water ice V and water ice VI. The square 
(coloured red in the online version) marks the critical point for Methane. 
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Fig. 3. — The phase diagram of methane filled ice-Ih water ice crystal. The diamond is a known point of 
stability ( Hirai et al.||2003" |. The four solid curves (blue in the on-line version) are the stability boundaries 
for the filled water ice, assuming for it a thermal expansivity equal to that for water ice VII and two, three 
and five times larger than that for water ice VII. A Lower thermal expansivity for filled water ice increases 
its stability, i.e. shifts the curve to the right to higher temperatures. The dash-dot curve is the melting curve 
for molecular water ice ( Lin et alT||2004 1 and the dashed curve is the melting curve for the super- ionic phase 
of water ( Goncharov et al!]|2009 1 . The square and circle data points represent the conditions at the base of 
the H/He atmosphere surrounding Neptune and Uranus, respectively (Redmer et al.|[20TT). 
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Fig. 4. — Dashed curve (red in the on-line version) is the enthalpy of fusion as a function of pressure for 



pure water (Goncharov et al. 20091. The branching in the melting curve, due to the super- ionic phase, is 



here at 47 GPa ( Goncharov et al.||2005 ) . Dashed-dotted curve (blue in the on-line version) is a hypothetical 
enthalpy of fusion, where, the branching in the melting curve occurs at 27 GPa. 
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Fig. 5. — The solid curve (blue) is the phase diagram for filled-ice Ih assuming for it a thermal expansivity 
twice as much as that of water ice VII and a branching in the melting curve at 47 GPa. The thick dashed 
curve (blue) is the variation in the phase diagram if the branching in the melting curve is shifted to 27 GPa. 
The dashed-dotted thick curves (green) confining the thick (blue) curve from both sides represent the change 
in the phase diagram if the enthalpy of fusion is varied by ±50% globally. 
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Fig. 6. — The probability of a methane molecule to occupy the large [yicHi 
of the structure I clathrate. 



and the small (y\cHi) cage 



